package cn.edu.neu.lab603

import java.util

import org.apache.commons.math3.stat.descriptive.DescriptiveStatistics
import org.apache.commons.math3.stat.regression.{OLSMultipleLinearRegression, SimpleRegression}

import scala.util.control.Breaks._

/** 一些常用的数学函数库。
  *
  * Created by yfwz100 on 2016/11/10.
  */
package object math {

  /**
    * Gets the median from a list of numbers.
    *
    * @param list the list of numbers
    * @return the median
    */
  def median(list: Array[Double]): Double = getStatistics(list).getPercentile(50)

  /**
    * Gets the Median absolute deviation (MAD) from a array of numbers.
    *
    * @param data the array of numbers
    * @return the mad
    */
  def mad(data: Array[Double]): Double = {
    var mad: Double = 0
    if (data.length > 0) {
      val median_ = median(data)
      val deviationSum: Array[Double] = data.map(d => Math.abs(median_ - d))
      mad = median(deviationSum)
    }
    mad
  }

  /**
    * Returns an object to compute descriptive statistics for an array of numbers.
    *
    * @param list the array of numbers. Must not be null.
    * @return descriptive statistics for the array of numbers.
    */
  def getStatistics(list: Array[Double]): DescriptiveStatistics = {
    // Get a DescriptiveStatistics instance
    new DescriptiveStatistics(list)
  }

  /**
    * Gets the Interquartile Range (IQR) from an array of numbers.
    *
    * @param data the array of numbers
    * @return the IQR
    */
  def iqr(data: Array[Double]): Double = {
    val s = data.sorted
    val q1 = (0.25 * (data.length + 1) - 1).toInt
    val q3 = (0.75 * (data.length + 1) - 1).toInt
    s(q3) - s(q1)
  }

  /**
    * Gets the length of the shortest row in a given matrix
    *
    * @param data the data matrix
    * @return the length of the shortest row int he matrix
    */
  def countShortestRow(data: Array[Array[Double]]): Int = {
    var minLength: Int = 0
    for (row <- data) {
      if (row.length < minLength) minLength = row.length
    }
    minLength
  }

  /**
    * Trims zeros at the end of an array.
    *
    * @param data the data array
    * @return the trimmed array
    */
  def trimZeroTail(data: Array[Double]): Array[Double] = util.Arrays.copyOfRange(data, 0, countNonZeroBeginning(data))

  /**
    * Counts the number of values different of zero at the beginning of
    * an array.
    *
    * @param data the array of numbers
    * @return the number of values different of zero at the beginning of the array
    */
  def countNonZeroBeginning(data: Array[Double]): Int = {
    var i: Int = data.length - 1
    breakable {
      while (i >= 0)
        if (data({
          i -= 1
          i + 1
        }) != 0) break
    }
    i + 2
  }

  /**
    * Gets the Local Regression (Loess) parameter estimates.
    *
    * @param y the y array
    * @return the Loess parameter estimates
    */
  def getLoessParameterEstimates(y: Array[Double]): Array[Double] = {
    val n: Int = y.length
    val x: Array[Double] = new Array[Double](n)
    for (i <- 0 until n) {
      x(i) = i + 1
    }
    createWeightedLinearRegression(x, y, getTricubeWeights(n)).regress.getParameterEstimates
  }

  def createWeightedLinearRegression(x: Array[Double], y: Array[Double], weights: Array[Double]): SimpleRegression = {
    val xW: Array[Double] = new Array[Double](x.length)
    val yW: Array[Double] = new Array[Double](y.length)
    var nubZeroWeights: Int = 0
    for (i <- weights.indices if weights(i) <= 0) {
      nubZeroWeights += 1
    }
    for (i <- x.indices) {
      if (nubZeroWeights >= 0.4 * weights.length) {
        // See: http://www.ncsu.edu/crsc/events/ugw07/Presentations/Crooks_Qiao/Crooks_Qiao_Alt_Presentation.pdf
        xW(i) = Math.sqrt(weights(i)) * x(i)
        yW(i) = Math.sqrt(weights(i)) * y(i)
      }
      else {
        xW(i) = x(i)
        yW(i) = y(i)
      }
    }
    createLinearRegression(xW, yW)
  }

  def createLinearRegression(x: Array[Double], y: Array[Double]): SimpleRegression = {
    val regression: SimpleRegression = new SimpleRegression
    for ((xi, yi) <- x zip y) {
      regression.addData(xi, yi)
    }
    regression
  }

  /**
    * Gets the tricube weights.
    *
    * @param n the number of weights
    * @return an array of tricube weights with n elements
    */
  def getTricubeWeights(n: Int): Array[Double] = {
    val weights: Array[Double] = new Array[Double](n)
    val top: Double = n - 1
    val spread: Double = top
    for (i <- 2 until n) {
      val k: Double = Math.pow(1 - Math.pow((top - i) / spread, 3), 3)
      if (k > 0)
        weights(i) = 1 / k
      else
        weights(i) = Double.MaxValue
    }
    weights(1) = weights(2)
    weights(0) = weights(1)
    weights
  }

  def createLinearRegression(x: Array[Array[Double]], y: Array[Double]): OLSMultipleLinearRegression = {
    val regression: OLSMultipleLinearRegression = new OLSMultipleLinearRegression
    regression.newSampleData(y, x)
    regression
  }

  /**
    * Gets the robust loess parameter estimates.
    *
    * @param y the y array
    * @return the robust loess parameter estimates
    */
  def getRobustLoessParameterEstimates(y: Array[Double]): Array[Double] = {
    val n: Int = y.length
    val x: Array[Double] = new Array[Double](n)
    for (i <- 0 until n) {
      x(i) = i + 1
    }
    val tricubeRegression: SimpleRegression = createWeightedLinearRegression(x, y, getTricubeWeights(n))
    val residuals: Array[Double] = new Array[Double](n)
    for (i <- 0 until n) {
      residuals(i) = y(i) - tricubeRegression.predict(x(i))
    }
    val tricubeBySquareRegression: SimpleRegression = createWeightedLinearRegression(x, y, getTricubeBisquareWeights(residuals))
    val estimates: Array[Double] = tricubeBySquareRegression.regress.getParameterEstimates
    if (java.lang.Double.isNaN(estimates(0)) || java.lang.Double.isNaN(estimates(1))) return tricubeRegression.regress.getParameterEstimates
    estimates
  }

  /**
    * Gets the tricube bisquare weights.
    *
    * @param residuals the residuals array
    * @return the tricube bisquare weights
    */
  def getTricubeBisquareWeights(residuals: Array[Double]): Array[Double] = {
    val n: Int = residuals.length
    val weights: Array[Double] = getTricubeWeights(n)
    val weights2: Array[Double] = new Array[Double](n)
    val s6: Double = median(residuals.map(Math.abs)) * 6
    for (i <- 2 until n) {
      val k: Double = Math.pow(1 - Math.pow(residuals(i) / s6, 2), 2)
      if (k > 0)
        weights2(i) = (1 / k) * weights(i)
      else
        weights2(i) = Double.MaxValue
    }
    weights2(1) = weights2(2)
    weights2(0) = weights2(1)
    weights2
  }

}
